* GECM DGP

// Basics
clear all
set matsize 11000
 
parallel setclusters 8

capture program drop ecm
program define ecm, rclass

// Monte Carlo Basics
version 12.1
syntax [, obs(integer 1) phi_u(real 0.5) beta_x(real 0.5) ]

drop _all

*args obs phi_u beta_x

set obs `obs'

// TS Basics
egen time = fill(0 1 2)
tsset time

// Data Generating Process
gen x = rnormal() if time == 0
replace x = l.x + rnormal() if time > 0

scalar gamma1 = 1.5
gen y = rnormal()
replace y = `phi_u'*l.y + `beta_x'*d.x + gamma1*l.x + rnormal() if time > 0

dfuller y
return scalar stat_y = r(p) < 0.05 

// Model estimation

* GECM
reg D.y L.y D.x L.x
return scalar phi_gecm = _b[L1.y]
local gecm_true_phi = `phi_u' - 1
scalar phi_gecm_test = (_b[L1.y] - `gecm_true_phi') / _se[L1.y]
return scalar phi_gecm_t1 = 2*ttail(e(df_r),abs(phi_gecm_test)) < 0.05
return scalar b1_gecm  = _b[D1.x]
scalar b1_gecm_test = (_b[D1.x] - `beta_x') / _se[D1.x]
return scalar b1_gecm_t1 = 2*ttail(e(df_r),abs(b1_gecm_test)) < 0.05
return scalar b2_gecm  = _b[L1.x]
scalar b2_gecm_test = (_b[L1.x] - (gamma1)) / _se[L1.x]
return scalar b2_gecm_t1 = 2*ttail(e(df_r),abs(b2_gecm_test)) < 0.05
return scalar lrm_gecm = -_b[L1.x] / _b[L1.y]
nlcom (-_b[L1.x] / _b[L1.y])* (1e+3)
matrix define var_lrm_gecm = r(V)
scalar var_lrm_gecm_se  = sqrt(var_lrm_gecm[1,1]) / (1e+3)
scalar lrm_gecm_test = ((-_b[L1.x] / _b[L1.y]) - ((-gamma1)/`gecm_true_phi')) /var_lrm_gecm_se
return scalar lrm_gecm_t1 = 2*ttail(e(df_r),abs(lrm_gecm_test)) < 0.05

scalar phi_gecm_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi_gecm_power = 2*ttail(e(df_r),abs(phi_gecm_test_power)) < 0.05
scalar b1_gecm_test_power = (_b[D1.x] - 0) / _se[D1.x]
return scalar b1_gecm_power = 2*ttail(e(df_r),abs(b1_gecm_test_power)) < 0.05
scalar b2_gecm_test_power = (_b[L1.x] - 0) / _se[L1.x]
return scalar b2_gecm_power = 2*ttail(e(df_r),abs(b2_gecm_test_power)) < 0.05
scalar lrm_gecm_test_power = ((-_b[L1.x] / _b[L1.y]) - 0) /var_lrm_gecm_se
return scalar lrm_gecm_power = 2*ttail(e(df_r),abs(lrm_gecm_test_power)) < 0.05


* Estimate ADL
reg y L.y x L.x
return scalar phi_adl = _b[L1.y]
scalar phi_adl_test = (_b[L1.y] - `phi_u') / _se[L1.y]
return scalar phi_adl_t1 = 2*ttail(e(df_r),abs(phi_adl_test)) < 0.05
return scalar b1_adl = _b[x]
scalar b1_adl_test = (_b[x] - `beta_x') / _se[x]
return scalar b1_adl_t1 = 2*ttail(e(df_r),abs(b1_adl_test)) < 0.05
return scalar b2_adl = _b[L1.x]
scalar b2_adl_test = (_b[L1.x] - (gamma1-`beta_x')) / _se[L1.x]
return scalar b2_adl_t1 = 2*ttail(e(df_r),abs(b2_adl_test)) < 0.05
return scalar lrm_adl = (_b[x] + _b[L1.x]) /(1-_b[L1.y])
nlcom ((_b[x] + _b[L1.x]) /(1-_b[L1.y])) * (1e+3)
matrix define var_lrm_adl = r(V)
scalar var_lrm_adl_se  = sqrt(var_lrm_adl[1,1]) / (1e+3)
scalar lrm_adl_test = (((_b[x] + _b[L1.x]) /(1-_b[L1.y])) - ((`beta_x' + (gamma1-`beta_x'))/ (1-`phi_u'))) /var_lrm_adl_se
return scalar lrm_adl_t1 = 2*ttail(e(df_r),abs(lrm_adl_test)) < 0.05

scalar phi_adl_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi_adl_power = 2*ttail(e(df_r),abs(phi_adl_test_power)) < 0.05
scalar b1_adl_test_power = (_b[x] - 0) / _se[x]
return scalar b1_adl_power = 2*ttail(e(df_r),abs(b1_adl_test_power)) < 0.05
scalar b2_adl_test_power = (_b[L1.x] - 0) / _se[L1.x]
return scalar b2_adl_power = 2*ttail(e(df_r),abs(b2_adl_test_power)) < 0.05
scalar lrm_adl_test_power = (((_b[x] + _b[L1.x]) /(1-_b[L1.y])) - 0) /var_lrm_adl_se
return scalar lrm_adl_power = 2*ttail(e(df_r),abs(lrm_adl_test_power)) < 0.05

* Estimate ADL (2,2) 
reg y L.y L2.y x L.x L2.x
return scalar phi1_adl22 = _b[L1.y]
scalar phi1_adl22_test = (_b[L1.y] - `phi_u') / _se[L1.y]
return scalar phi1_adl22_t1 = 2*ttail(e(df_r),abs(phi1_adl22_test)) < 0.05
return scalar phi2_adl22 = _b[L2.y]
scalar phi2_adl22_test = (_b[L2.y] - 0) / _se[L2.y]
return scalar phi2_adl22_t1 = 2*ttail(e(df_r),abs(phi2_adl22_test)) < 0.05
return scalar b1_adl22 = _b[x]
scalar b1_adl22_test = (_b[x] - `beta_x') / _se[x]
return scalar b1_adl22_t1 = 2*ttail(e(df_r),abs(b1_adl22_test)) < 0.05
return scalar b2_adl22 = _b[L1.x]
scalar b2_adl22_test = (_b[L1.x] - (gamma1-`beta_x')) / _se[L1.x]
return scalar b2_adl22_t1 = 2*ttail(e(df_r),abs(b2_adl22_test)) < 0.05
return scalar b3_adl22 = _b[L2.x]
scalar b3_adl22_test = (_b[L2.x] - 0) / _se[L2.x]
return scalar b3_adl22_t1 = 2*ttail(e(df_r),abs(b3_adl22_test)) < 0.05
return scalar lrm_adl22 = (_b[x] + _b[L1.x] + _b[L2.x]) /(1-_b[L1.y]-_b[L2.y])
nlcom ((_b[x] + _b[L1.x] + _b[L2.x]) /(1-_b[L1.y]-_b[L2.y]))*(1e+3)
matrix define var_lrm_adl22 = r(V)
scalar var_lrm_adl22_se  = sqrt(var_lrm_adl22[1,1]) / (1e+3)
scalar lrm_adl22_test = (((_b[x] + _b[L1.x] + _b[L2.x]) /(1-_b[L1.y]-_b[L2.y])) - ((`beta_x' + (gamma1-`beta_x'))/ (1-`phi_u'))) /var_lrm_adl22_se
return scalar lrm_adl22_t1 = 2*ttail(e(df_r),abs(lrm_adl22_test)) < 0.05

scalar phi1_adl22_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi1_adl22_power = 2*ttail(e(df_r),abs(phi1_adl22_test_power)) < 0.05
scalar b1_adl22_test_power = (_b[x] - 0) / _se[x]
return scalar b1_adl22_power = 2*ttail(e(df_r),abs(b1_adl22_test_power)) < 0.05
scalar b2_adl22_test_power = (_b[L1.x] - 0) / _se[L1.x]
return scalar b2_adl22_power = 2*ttail(e(df_r),abs(b2_adl22_test_power)) < 0.05
scalar lrm_adl22_test_power = (((_b[x] + _b[L1.x] + _b[L2.x]) /(1-_b[L1.y]-_b[L2.y])) - 0) /var_lrm_adl22_se
return scalar lrm_adl22_power = 2*ttail(e(df_r),abs(lrm_adl22_test_power)) < 0.05

* Estimate ADL (3,3) 
reg y L.y L2.y L3.y x L.x L2.x L3.x
return scalar phi1_adl33 = _b[L1.y]
scalar phi1_adl33_test = (_b[L1.y] - `phi_u') / _se[L1.y]
return scalar phi1_adl33_t1 = 2*ttail(e(df_r),abs(phi1_adl33_test)) < 0.05
return scalar phi2_adl33 = _b[L2.y]
scalar phi2_adl33_test = (_b[L2.y] - 0) / _se[L2.y]
return scalar phi2_adl33_t1 = 2*ttail(e(df_r),abs(phi2_adl33_test)) < 0.05
return scalar phi3_adl33 = _b[L3.y]
scalar phi3_adl33_test = (_b[L3.y] - 0) / _se[L3.y]
return scalar phi3_adl33_t1 = 2*ttail(e(df_r),abs(phi3_adl33_test)) < 0.05
return scalar b1_adl33 = _b[x]
scalar b1_adl33_test = (_b[x] - `beta_x') / _se[x]
return scalar b1_adl33_t1 = 2*ttail(e(df_r),abs(b1_adl33_test)) < 0.05
return scalar b2_adl33 = _b[L1.x]
scalar b2_adl33_test = (_b[L1.x] - (gamma1-`beta_x')) / _se[L1.x]
return scalar b2_adl33_t1 = 2*ttail(e(df_r),abs(b2_adl33_test)) < 0.05 
return scalar b3_adl33 = _b[L2.x]
scalar b3_adl33_test = (_b[L2.x] - 0) / _se[L2.x]
return scalar b3_adl33_t1 = 2*ttail(e(df_r),abs(b3_adl33_test)) < 0.05
return scalar b4_adl33 = _b[L3.x]
scalar b4_adl33_test = (_b[L3.x] - 0) / _se[L3.x]
return scalar b4_adl33_t1 = 2*ttail(e(df_r),abs(b4_adl33_test)) < 0.05
return scalar lrm_adl33 = (_b[x] + _b[L1.x] + _b[L2.x] + _b[L3.x]) /(1-_b[L1.y]-_b[L2.y]-_b[L3.y])
nlcom ((_b[x] + _b[L1.x] + _b[L2.x] + _b[L3.x]) /(1-_b[L1.y]-_b[L2.y]-_b[L3.y])) * (1e+3)
matrix define var_lrm_adl33 = r(V)
scalar var_lrm_adl33_se  = sqrt(var_lrm_adl33[1,1]) / (1e+3)
scalar lrm_adl33_test = (((_b[x] + _b[L1.x] + _b[L2.x] + _b[L3.x]) /(1-_b[L1.y]-_b[L2.y]-_b[L3.y])) - ((`beta_x' + (gamma1-`beta_x'))/ (1-`phi_u'))) /var_lrm_adl33_se
return scalar lrm_adl33_t1 = 2*ttail(e(df_r),abs(lrm_adl33_test)) < 0.05

scalar phi1_adl33_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi1_adl33_power = 2*ttail(e(df_r),abs(phi1_adl33_test_power)) < 0.05
scalar b1_adl33_test_power = (_b[x] - 0) / _se[x]
return scalar b1_adl33_power = 2*ttail(e(df_r),abs(b1_adl33_test_power)) < 0.05
scalar b2_adl33_test_power = (_b[L1.x] - 0) / _se[L1.x]
return scalar b2_adl33_power = 2*ttail(e(df_r),abs(b2_adl33_test_power)) < 0.05
scalar lrm_adl33_test_power = (((_b[x] + _b[L1.x] + _b[L2.x] + _b[L3.x]) /(1-_b[L1.y]-_b[L2.y]-_b[L3.y]))  - 0) /var_lrm_adl33_se
return scalar lrm_adl33_power = 2*ttail(e(df_r),abs(lrm_adl33_test_power)) < 0.05

* F-test, ADL(2,2)
reg y L.y L2.y x L.x L2.x
test L.y + L2.y = `phi_u'
return scalar ftest_y_adl22_t1 = r(p) < 0.05
reg y L.y L2.y x L.x L2.x
test L.y + L2.y = 0
return scalar ftest_y_adl22_power = r(p) < 0.05
reg y L.y L2.y x L.x L2.x
test x + L.x + L2.x = `beta_x' + gamma1-`beta_x'
return scalar ftest_x_adl22_t1 = r(p) < 0.05
reg y L.y L2.y x L.x L2.x
test x + L.x + L2.x = 0
return scalar ftest_x_adl22_power = r(p) < 0.05


* F-test, ADL(3,3)
reg y L.y L2.y L3.y x L.x L2.x L3.x
test L.y + L2.y + L3.y = `phi_u'
return scalar ftest_y_adl33_t1 = r(p) < 0.05
reg y L.y L2.y L3.y x L.x L2.x L3.x
test L.y + L2.y + L3.y = 0
return scalar ftest_y_adl33_power = r(p) < 0.05
reg y L.y L2.y L3.y x L.x L2.x L3.x
test x + L.x + L2.x + L3.x = `beta_x' + gamma1-`beta_x'
return scalar ftest_x_adl33_t1 = r(p) < 0.05
reg y L.y L2.y L3.y x L.x L2.x L3.x
test x + L.x + L2.x + L3.x = 0
return scalar ftest_x_adl33_power = r(p) < 0.05


end
